package edu.hawaii.nearestneighbor.util;

/**
 * Provides the functions for estimating gamma.<br>
 * From sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN
 * 0-521-43108-5) Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C)
 * 1988-1992 by Numerical Recipes Software. Permission is granted for internet users to make one
 * paper copy for their own personal use. Further reproduction, or any copying of machinereadable
 * files (including this one) to any server computer, is strictly prohibited. To order Numerical
 * Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America
 * only), or send email to directcustserv@cambridge.org (outside North America).
 * See: <a href="http://www.library.cornell.edu/nr/bookcpdf/c6-1.pdf">Link to source</a>.
 * 
 * @author Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley).
 * @author Hastings, C. 1955, Approximations for Digital Computers (Princeton: Princeton University
 *         Press).
 * @author Luke, Y.L. 1975, Mathematical Functions and Their Approximations (New York: Academic
 *         Press).
 * 
 */
public final class Gamma {

  /**
   * Private constructor.
   */
  private Gamma() {
    // empty
  }
  
  /**
   * Calculates the estimate to gamma(x).
   * @param x The parameter x in the gamma function.
   * @return An estimation to gamma(x).
   */
  public static double gamma(double x) {
    return Math.exp(logGamma(x));
  }

  /**
   * Calculates the log of gamma.
   * Adapted from C language.
   * @param xx The parameter x in the gamma function.
   * @return An estimation to ln of gamma(x).
   */
  public static double logGamma(double xx) {
    double x, y, tmp, ser;
    double[] cof = { 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155,
        0.1208650973866179e-2, -0.5395239384953e-5 };
    int j;
    x = xx;
    y = x;
    tmp = x + 5.5;
    tmp -= (x + 0.5) * Math.log(tmp);
    ser = 1.000000000190015;
    for (j = 0; j <= 5; j++) {
      ser += cof[j] / ++y;
    }
    return -tmp + Math.log(2.5066282746310005 * ser / x);
  }

}
